clear 
est clear

********************************************************************************
*	Programs
******************************************************************************** 

* program for the OLS regressions... 
capture program drop lm_regs
program define lm_regs, eclass

	*** Syntax 
	syntax varlist(numeric) [, y_bl(varlist numeric) ctrls(varlist numeric) treatment(varlist numeric) rounds(numlist) t_suff(string)]

	*** Unadjusted Regressions 
qui	{
	foreach y of local varlist {
		foreach r in `rounds' {
		*** Create the centered block variables 
			capture drop blocks*
			qui tab block, gen(blocks)
			qui foreach v of varlist blocks* {
					replace `v' = . 			if missing(`y') & round ==`r'
					su `v' 						if round == `r' & !missing(`y')
					replace `v' = `v' - r(mean) if round == `r' & !missing(`y')
				}
		*** Linear Model: treatment, Block FE, Cluster SEs at Club level, weights (tracking) 
		eststo `y'_`t_suff'_`r':		reg `y' b0.`treatment'##c.(blocks*) [pw=w] if round==`r', vce(cluster club)
			qui {
				*** add scalars
				* control mean
				qui su `y' [aw=w] if `treatment'==0 & round==`r' 
				estadd scalar c_mean = r(mean) 	
				estadd scalar c_sd = r(sd)
				* p-value 
				test 1.`treatment' == 0 
				estadd scalar pval = r(p) 
			}
		}
	}
}
	**** Adjusted Regressions 
qui {
	foreach y of local varlist {
		foreach r in `rounds' {
		*** Create the centered block variables 
			capture drop blocks*
			qui tab block, gen(blocks)
			qui foreach v of varlist blocks* {
					replace `v' = . 			if missing(`y')
					su `v' 						if round == `r' & !missing(`y')
					replace `v' = `v' - r(mean) if round == `r' & !missing(`y')
				}
		*** Linear Model: treatment, Block FE, Cluster SEs at Club level, weights (tracking), BL covars fully interacted
		eststo `y'_`t_suff'_`r'_adj:	reg `y' b0.`treatment'##c.(`ctrls' `y_bl' blocks*) [pw=w] if round==`r', vce(cluster club)
			qui {
				*** add scalars
				* control mean
				qui su `y' [aw=w] if `treatment'==0 & round==`r'   
				estadd scalar c_mean = r(mean) 	
				estadd scalar c_sd = r(sd)
				* p-value 
				test 1.`treatment' == 0 
				estadd scalar pval = r(p) 
			}
		}
	}
}
end


* program to extract q values for the outcomes specified BY round
capture program drop get_qvals
program define get_qvals, eclass

	*** Syntax 
	syntax varlist(numeric) [, rounds(numlist) adj(string) t_suff(string) method(string) iv(string)]

	*** if LATE change the suffix 
	if ("`iv'"=="yes") local t_suff = "iv"
	di "`t_suff'"
	
	*** Create vector of p-values
	* matrix to store them
	local varcount : word count `varlist'
	matrix pvals_tab = J(`varcount',1,.) 
	matrix colnames pvals_tab = pvals
	* loop over outcomes and rounds
	qui foreach r in `rounds' {
	    local i = 0
		foreach y of local varlist {
		    local ++i
			* restore estimates
			est restore `y'_`t_suff'_`r'`adj' 
			est replay
			* collect pvals in matrix 
			if ("`iv'"=="yes") {
				local beta = e(b)[1,"1.treat_ther"]
				local beta_se = sqrt(e(V)["1.treat_ther","1.treat_ther"])
				local beta_z = `beta'/`beta_se'
				local beta_p = 2*(1-normal(abs(`beta_z')))
				matrix pvals_tab[`i',1] = `beta_p'
				matrix list pvals_tab
			}
			else if ("`iv'"!="yes") {
				matrix list r(table)
				matrix pvals_tab[`i',1] = r(table)["pvalue","1.treat_`t_suff'"]
			}
			* save in a dataset 
			capture drop pvals`r'
			svmat pvals_tab, names(col)
			rename pvals pvals`r'
			
			*** get q-values  
			capture drop qvals`r'
			qqvalue pvals`r', method(`method') qvalue(qvals`r') 
		}
	}
	*** add to estimates 
	* loop over outcomes and rounds
	qui foreach r in `rounds' {
	    local i = 0
		foreach y of local varlist {
		    local ++i
			* restore estimates
			est restore `y'_`t_suff'_`r'`adj' 
			* add as a scalar
			local qv = qvals`r' in `i'
			di `qv'
			estadd scalar qval = `qv'
		}
	} 
end






********************************************************************************
*	I - Impacts of Therapy vs. Control
********************************************************************************

use "${data}/Analysis_wide.dta", clear
distinct block club

*** is it easier in long mode? I think it might be...
reshape long phq8_score ghq12_score rosb_score res_score loc_score ghq12_min phq8_min w, i(block cr_id club treat ) j(round)

*** drop +cash folks in longer follow ups 
drop if treat==2 & round>=2
tab round 

*** run regs for main outcomes  
lm_regs phq8_min phq8_score ghq12_min ghq12_score , ///
	treatment(treat_ther) t_suff("ther") y_bl(phq8_score_b)  ///
	ctrls(age_b everpregnant_b evermarried_b PPI_score_b) rounds(1 2 3)

*** run regs for secondary
lm_regs rosb_score res_score loc_score, ///
	treatment(treat_ther) t_suff("ther") y_bl(phq8_score_b) ///
	ctrls(age_b everpregnant_b evermarried_b PPI_score_b) rounds(1 3) 

*** compute q-values 
* "simes" is the Benjamini Hochberg (1995) step-up procedure
get_qvals phq8_min ghq12_min , rounds(1 2 3) adj(" ") t_suff("ther") method("simes")
get_qvals phq8_min ghq12_min , rounds(1 2 3) adj("_adj") t_suff("ther") method("simes")

get_qvals phq8_score ghq12_score rosb_score res_score loc_score , rounds(1 3) adj(" ") t_suff("ther") method("simes")
get_qvals phq8_score ghq12_score rosb_score res_score loc_score , rounds(1 3) adj("_adj") t_suff("ther") method("simes")

get_qvals phq8_score ghq12_score  , rounds(2) adj(" ") t_suff("ther") method("simes")
get_qvals phq8_score ghq12_score  , rounds(2) adj("_adj") t_suff("ther") method("simes")



 	
 
*** unadjusted table
* RR 
#delimit ;
esttab 	phq8_min_ther_1 ghq12_min_ther_1 
		phq8_score_ther_1 ghq12_score_ther_1 	 
		rosb_score_ther_1 res_score_ther_1 loc_score_ther_1
	using "${tables}/results-mh-both1.tex", replace  booktabs
	label nostar nonote nomtitle nonumber frag nolines nogaps
	b(3) se(3) keep(1.treat_ther) 
	stats(c_mean c_sd N pval qval, fmt(3 3 0 3 3)  
	label("Control mean" "Control SD" "Observations" "H0: IPT-G=0 p-values" "\hspace{20pt}FDR adj. q-values")) 
	; 
#delimit cr 
* Midline 
#delimit ;
esttab 	phq8_min_ther_2   ghq12_min_ther_2
		phq8_score_ther_2 ghq12_score_ther_2 	
	using "${tables}/results-mh-both2.tex", replace  booktabs
	label nostar nonote nomtitle nonumber frag nolines nogaps
	b(3) se(3) keep(1.treat_ther) 
	stats(c_mean c_sd N pval qval, fmt(3 3 0 3 3)  
	label("Control mean" "Control SD" "Observations" "H0: IPT-G=0 p-values" "\hspace{20pt}FDR adj. q-values")) 
	; 
#delimit cr 
* Endline 
#delimit ;
esttab 	phq8_min_ther_3  ghq12_min_ther_3  	
		phq8_score_ther_3 ghq12_score_ther_3 	 
		rosb_score_ther_3 res_score_ther_3 loc_score_ther_3
	using "${tables}/results-mh-both3.tex", replace  booktabs
	label nostar nonote nomtitle nonumber frag nolines nogaps
	b(3) se(3) keep(1.treat_ther) 
	stats(c_mean c_sd N pval qval, fmt(3 3 0 3 3)  
	label("Control mean" "Control SD" "Observations" "H0: IPT-G=0 p-values" "\hspace{20pt}FDR adj. q-values")) 
	; 
#delimit cr 

*** adjusted table
* RR 
#delimit ;
esttab 	phq8_min_ther_1_adj ghq12_min_ther_1_adj 
		phq8_score_ther_1_adj ghq12_score_ther_1_adj 	 
		rosb_score_ther_1_adj res_score_ther_1_adj loc_score_ther_1_adj
	using "${tables}/results-mh-both1-adj.tex", replace  booktabs
	label nostar nonote nomtitle nonumber frag nolines nogaps
	b(3) se(3) keep(1.treat_ther) 
	stats(c_mean c_sd N pval qval, fmt(3 3 0 3 3)  
	label("Control mean" "Control SD" "Observations" "H0: IPT-G=0 p-values" "\hspace{20pt}FDR adj. q-values")) 
	; 
#delimit cr 
* Midline 
#delimit ;
esttab 	phq8_min_ther_2_adj   ghq12_min_ther_2_adj
		phq8_score_ther_2_adj ghq12_score_ther_2_adj 	
	using "${tables}/results-mh-both2-adj.tex", replace  booktabs
	label nostar nonote nomtitle nonumber frag nolines nogaps
	b(3) se(3) keep(1.treat_ther) 
	stats(c_mean c_sd N pval qval, fmt(3 3 0 3 3)  
	label("Control mean" "Control SD" "Observations" "H0: IPT-G=0 p-values" "\hspace{20pt}FDR adj. q-values")) 
	; 
#delimit cr 
* Endline 
#delimit ;
esttab 	phq8_min_ther_3_adj ghq12_min_ther_3_adj 	
		phq8_score_ther_3_adj ghq12_score_ther_3_adj 	 
		rosb_score_ther_3_adj res_score_ther_3_adj loc_score_ther_3_adj
	using "${tables}/results-mh-both3-adj.tex", replace  booktabs
	label nostar nonote nomtitle nonumber frag nolines nogaps
	b(3) se(3) keep(1.treat_ther) 
	stats(c_mean c_sd N pval qval, fmt(3 3 0 3 3)  
	label("Control mean" "Control SD" "Observations" "H0: IPT-G=0 p-values" "\hspace{20pt}FDR adj. q-values")) 
	; 
#delimit cr 

*/








 














********************************************************************************
*	III - Impacts of Cash + Therapy over just Therapy
********************************************************************************

est clear

use "${data}/Analysis_wide.dta", clear
distinct block club

*** is it easier in long mode? I think it might be...
reshape long phq8_score ghq12_score rosb_score res_score loc_score ghq12_min phq8_min w, i(block cr_id club treat ) j(round)

*** drop +cash folks in longer follow ups 
drop if treat==0 
drop if round==1
tab round 

*** value label of cash arm 
label define treat_cash 0 "IPT-G" 1 "IPT-G+", modify

*** run regs for main outcomes 
lm_regs phq8_min phq8_score ghq12_min ghq12_score , rounds(2 3) ///
	treatment(treat_cash) t_suff("cash") y_bl(phq8_score_b) ///
	ctrls(age_b everpregnant_b evermarried_b PPI_score_b) 
	
lm_regs rosb_score res_score loc_score , rounds(3) ///
	treatment(treat_cash) t_suff("cash") y_bl(phq8_score_b) ///
	ctrls(age_b everpregnant_b evermarried_b PPI_score_b)



*** compute q-values 
* "simes" is the Benjamini Hochberg (1995) step-up procedure
get_qvals phq8_min ghq12_min , rounds(2 3) adj(" ") t_suff("cash") method("simes")
get_qvals phq8_min ghq12_min , rounds(2 3) adj("_adj") t_suff("cash") method("simes")

get_qvals phq8_score ghq12_score rosb_score res_score loc_score , rounds(3) adj(" ") t_suff("cash") method("simes")
get_qvals phq8_score ghq12_score rosb_score res_score loc_score , rounds(3) adj("_adj") t_suff("cash") method("simes")

get_qvals phq8_score ghq12_score  , rounds(2) adj(" ") t_suff("cash") method("simes")
get_qvals phq8_score ghq12_score  , rounds(2) adj("_adj") t_suff("cash") method("simes")



*** unadjusted table
* Midline
#delimit ;
esttab 	phq8_min_cash_2   ghq12_min_cash_2
		phq8_score_cash_2 ghq12_score_cash_2 	
	using "${tables}/results-mh-both2-cash.tex", replace  booktabs
	label nostar nonote nomtitle nonumber frag nolines nogaps
	b(3) se(3) keep(1.treat_cash) 
	stats(c_mean c_sd N pval qval, fmt(3 3 0 3 3)  
	label("IPT-G mean" "IPT-G SD" "Observations" "H0: IPT-G+=0 p-values" "\hspace{20pt}FDR adj. q-values")) 
	; 
#delimit cr 
* Endline
#delimit ;
esttab 	phq8_min_cash_3 ghq12_min_cash_3  	
		phq8_score_cash_3 ghq12_score_cash_3 	 
		rosb_score_cash_3 res_score_cash_3 loc_score_cash_3
	using "${tables}/results-mh-both3-cash.tex", replace  booktabs
	label nostar nonote nomtitle nonumber frag nolines nogaps
	b(3) se(3) keep(1.treat_cash) 
	stats(c_mean c_sd N pval qval, fmt(3 3 0 3 3)  
	label("IPT-G mean" "IPT-G SD" "Observations" "H0: IPT-G+=0 p-values" "\hspace{20pt}FDR adj. q-values")) 
	; 
#delimit cr 

*** adjusted table
* Midline
#delimit ;
esttab 	phq8_min_cash_2_adj   ghq12_min_cash_2_adj
		phq8_score_cash_2_adj ghq12_score_cash_2_adj 	
	using "${tables}/results-mh-both2-cash-adj.tex", replace  booktabs
	label nostar nonote nomtitle nonumber frag nolines nogaps
	b(3) se(3) keep(1.treat_cash) 
	stats(c_mean c_sd N pval qval, fmt(3 3 0 3 3)  
	label("IPT-G mean" "IPT-G SD" "Observations" "H0: IPT-G+=0 p-values" "\hspace{20pt}FDR adj. q-values")) 
	; 
#delimit cr 
* Endline
#delimit ;
esttab 	phq8_min_cash_3_adj  ghq12_min_cash_3_adj 	
		phq8_score_cash_3_adj  ghq12_score_cash_3_adj 	 
		rosb_score_cash_3_adj res_score_cash_3_adj loc_score_cash_3_adj
	using "${tables}/results-mh-both3-cash-adj.tex", replace  booktabs
	label nostar nonote nomtitle nonumber frag nolines nogaps
	b(3) se(3) keep(1.treat_cash) 
	stats(c_mean c_sd N pval qval, fmt(3 3 0 3 3)  
	label("IPT-G mean" "IPT-G SD" "Observations" "H0: IPT-G+=0 p-values" "\hspace{20pt}FDR adj. q-values")) 
	; 
#delimit cr 

